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■i-T 1 INTRODUCTION 



ABSTRACT 

We present measurements of the number density of voids in the dark matter distribu- 
tion from a series of N-body simulations of a ACDM cosmology. We define voids as 
spherical regions of p v = 0.2p m around density minima in order to relate our results to 
the predicted abundances using the excursion set formalism. Using a linear underden- 
sity of S v = —2.7, from a spherical evolution model, we find that a volume conserving 
model, which does not conserve number density in the mapping from the linear to 
nonlinear regime, matches the measured abundance to within 16% for a range of void 
radii 1 < r(/i _1 Mpc) < 15. This model fixes the volume fraction of the universe which 
is in voids and assumes that voids of a similar size merge as they expand by a factor 
of 1.7 to achieve a nonlinear den sity of p v = 0.2p m today. We find that the model of 
iSheth fc van de Wevgaertl (|2004[ ) for the number density of voids greatly overpredicts 
the abundances over the same range of scales. We find that the volume conserving 
model works well at matching the number density of voids measured from the simula- 
tions at higher redshifts, z = 0.5 and 1, as well as correctly predicting the abundances 
to within 25% in a simulation of a matter dominated fl m — 1 universe. We examine 
the abundance of voids in the halo distribution and find fewer small, r < 10/i _1 Mpc, 
voids and many more large, r > 10/i -1 Mpc, voids compared to the dark matter. These 
results indicate that voids identified in the halo or galaxy distribution are related to 
the underlying void distribution in the dark matter in a complicated way which merits 
further study if voids are to be used as a precision probe of cosmology. 

Key words: Methods: iV-body simulations - Cosmology: theory - large-scale struc- 
ture of the Universe 



Galaxy redshift surveys allow us to study and map out 
the large scale structure of our Universe revealing a hier- 
archical mass distribution with substructure over a wide 
range of scales. The main components of the galaxy 
distribution are arranged in a r emarkable 'cosmic web' 
(jBond. Kofman fc Pogosvanl 1 19961 ) made up of clusters of 
galaxies connected by filaments with large empty voids 
which occupy most of the volume. Only recently have 
systematic studies using voids as precision probes of the 
growth of structure been possible due to the increased depth 
and volume of current galaxy surveys (IColless et al.ll200ll ; 
lYork. fc SPSS Collaboration! l200d : lAbazaiian et all (20091 ). 
In this paper we study the distribution of underdense void 
regions in the dark matter and halo distributions using N- 
body simulations. We focus on the excursion set method 
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which gives an analytical prescription for the number den- 
sity of voids which we compare with measurements from 
simulations. 

Voids are a common feature in galaxy surveys 
with one of the most well known discoveries being the 
void in Bo otes which has a d iameter of approximately 
50/i _1 Mpc (jKirshner et al.l 119811 ). Since then several sur- 
veys such as the Cente r for Astrophysics Redshift Survey 
' Geller fc Huchrall 19891), the Southern Sky Redshift Survey 



Maurogordato. Schaeffer fc da Costal 1992T) and the d eeper 



Las Campanas Redshift Survey ( Shectman et al.ll 19961 ) have 
identified voids in the distribution of galaxies and clusters 
confirming that they are the dominant a nd volume filling 
com ponent of our Univer se. Most recently |Pan et al.l (|2012l ) 
and ISutter et ail l|2012bl ) both used the Sloan Digital Sky 
Survey Data Rele ase 7 (SPSS DR7 ) JAbazaiian et alj|200a ) 
to identify voids. |Pan et al.l (|2012T ) found 1054 statistically 
significant voids with radii r > 10/i _1 Mpc with an absolute 
magnitude cut of M r < —20.09. They argue that voids of 
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effective radius r c g ~ 20ft~ 1 Mpc dominate the void volume 
with the largest void in the sample having r ~ 30/i _1 Mpc. 
ISutter et al.l (|2012bl ) constructed the first public void cat- 
alog using the full extent of the SDSS DR7 spectroscopic 
survey which included the LRGs and found large voids in 
the sample of r ~ 50 — 60/i _1 Mpc in radius. 

Early numerical and theoret ic al work on the ev o lution 
of vo ids bvlRegos fc Gelled (|l99ll ). lBlumenthal et al.l (fl993 ) 



and iDubinski et al.l ( 19931 ) focused on the expansion of 
initial linear undensities up to the moment of shell cross- 
ing, which is used to define a characteristic time in the 
formation of voids. They considered spherical voids in a 
fi m = 1 universe and found that shell crossing occurs at a 
linear underdensity of —2.7 at which point the comoving 
size of the void has increased by a factor of 1.7 (see also 
van de Weygaert fc van Kampenlll993l; iFriedmann fc Piranl 



20011 ) . Many studies since then have focused on analysing 



the dynamics and statistical properties of voids such as 
the void probability function (VPF), the probabilit y that 
a ran domly placed sphere will contain no objects (jWhitel 
1 19791 ). the void filling factor, th e void number density 
profi l es (see e.g van de Weygaert 199ll; 



and void density 



cy pron . 
Mathis fc White! |2002|; iBenson et al. . 
20051; Shandarin et al.l 12006 ; Betancort-Riio et al 



Einasto et al 



Kreckel et al 



20031; 



Colberg et al.l 



201 ll; Ivan de Weygaert fc Plat 
l2012l ; lAragon-Calvo fc Szalavll201 ' 



2009; 



2011: 



Recent studies have looked at stacking voids in order to 
increase the statistical sign i ficanc e of weak lens i ng sig nals 
Jfflguchi, Qguri fc Hamanal |2012|; iKrause et al] l2013T )(see 
also lAmendola. Frieman fc Waga _ 19991) . the Integrated 
Sach - Wolfe effect ( Ilic. Langer fc Douspisl 120131 : ICai et al.l 
|2013J) or to extract cosmological parameters by modelling 
the distortions in redshift space (jLavaux fc WandeltJ 120121 : 
Sutter et al. 2012a; Bos ct al. 2012) or as a test of modified 
gravity ( Clampitt. Cai fc Lill2013l ). The precision of these 
tests relies on many factors, for example, given a survey 
or numerical simulation of a certain size, how robustly can 
we measure statistics for a void of a given size; how accu- 
rately can we predict the number density of voids in the 
galaxy/dark matter halo distribution and how well do voids 
in the galaxy/halo distribution trace voids in the dark mat- 
ter. In this paper we address these three issues. In our dis- 
cussion of a robust void finder we do not compare with all 
the other algorithms which have been used in previous stud- 
ies - our choice of void finder is motivated by the excursion 
set formalism for the abundance of voids which we aim to 
test. 

In analysing both galaxy surveys and numerical 
simulations a wide variety of void finding algorithms 
have been us e d to define underdense regions as voids. 
IColberg et al.l (|2008l ) carried out the first systematic re- 
view of 13 different void finders, identifying only two ar- 
eas of agreement amongst the different algorithms; that 
voids are very underdense (p ~ 0.05/9 m ) at their cen- 
tres and that voids have very steep edges. The void 
finding methods include the construction of proto-voids 
around local minima in the smoothed density field, af- 
ter separating the galaxy sample into ' wall' and 'void' 
galax ies (see e.g. |E1-Ad fc Piranl ll997J ; iHovle fc Vogelevl 
120021 ); merging proto-voids which results in non-spherical 
voids (|Colberg et al] 120051): the wat ershed algorithm 
(|Platen. van de Weygaert fc Jonesl 120071 ') which uses the 



DTF E method (|Schaap| 120071 ; ICautun fc van de Weygaertl 
120111 ). The watershed void finder identifies minima in the 
density field and construct voids by flooding basins until the 
'landscape' resembles a segmented plane where the edges 
of each segment outline a void region. A similar algorithm, 
which we make use of in this paper, is the ZOBOV iJNevrinckl 
120081 ) void finde r which uses the Voron oi tessellation field 
method (see e.g. Ivan de Weygaertl 120071 ) to partition parti- 
cles into zones, which are then joined together near density 
minima, into voids. 



There have been many studies of the excursion set 
method to predict th e abundance o f dark matter halos in 
the Universe (see e.g. IZentnerll2007l . for a review). In com- 
parison fewer studies have focused on testing the excursion 
set predictions for underdensities in the dark matter dis- 
tribution and we briefly outline some of these works here. 
In app lying this method to voids. ISheth fc van de Weygaertl 
((200J) presented a model for the abundance of voids in 
the dark matter which included the influence of the larger 
scale environment on the formation of a void. Their model 
takes into account two effects, firstly, a void of a given 
size may be embedded in another underdense region which 
is on a larger scale, the 'void-in-void' scenario, and sec- 
ondly, a void of a given size could be embedded in an 
overd ense region on a larger sca le, the 'void-in-cloud' sce- 
nario. iFurlanetto fc Piranl 1)20061 ) analysed how the barrier 
for shell crossing of a void in the galaxy distribution would 
differ from the linear theory barrier for dark matter, find- 
ing that voids selected from catalogs of luminous galax- 
ies should b e larger than those selected f r om faint galax- 
ies (se e also iD'Aloisio fc Furlanettol 120071 ). iD'Amico et al] 
(|2011f ) consider using voids as a probe of primordial non- 
Gaussianity and calculate the abundance of voids using the 
e xcursion set formalism and the t wo barrier presc r iption 
of lSheth fc van de Weygaertl (|2004l ) . IShandarin ei~aH (|2006l ) 
define voids as isolated regions of the low-density excur- 
sion set specified by density thresholds and measured the 
abundance and morphology of voids using N-body simu- 
lati ons. In this paper we wish to test the predictions of 
the lSheth fc van de Weygaertl (|2004l ) excursion set model by 
comparing them to measurements of void abundances from 
N-body simulati ons. Our definition of a void is similar to 
that adopted by IShandarin et al.l (|2006h as we use a strict 
density threshold to define the void edge (although we do 
not use isodensity contours). To our knowledge, this is the 
first time that this model has been directly compared with 
numerical simulations. 



This paper is organised as follows. In Section [2] we dis- 
cuss the excursion set method as it applies to dark matter 
halos and voids. Appendices [X] and [B] review the salient fea- 
tures of the spherical evolution model that connects the two. 
In Section [3] we detail our void finder and the N-body sim- 
ulations that were carried out. In Section U we present the 
main results of this paper on the number density of voids in 
three different cosmological models at 2 = and we show 
how this abundance changes with redshift. We also present 
the measured number density of voids in the halo distribu- 
tion. In Section [5] we present our conclusions. 



© 2013 RAS, MNRAS 000.[TIH6l 



The abundance of voids and the excursion set formalism 



2 EXCURSION SET FORMALISM FOR VOID 
ABUNDANCE 

In this section we begin with a brief review of the excur- 
sion set formalism in Section 12.11 It is well known that in 
combination with spherical collapse this approach provides 
insight into many aspects of halo formation and can be used 
to predict dark m atter halo abundances and clustering (see 
e.g |Zentnerll2007l . for a review of the subject). The analogous 
spherical expansion model can li kewise be used to make ex- 
cursi on set predictions for voids (jSheth fc van de Weygaertl 
[200J). We review this extension in Section [2~2l and show that 
it requires modifications on physical grounds. We propose a 
simple modification based on volume fraction conservation 
in Section [2731 



2.1 Excursion set formalism 

The excursion set formalism at its heart relies on knowledge 
of the statistical properties of the linear density field. In 
Fourier space, the linear density fluctuation field smoothed 
on a scale R is given by 



5(2, R) = 



d k S(k)W(k,R)e-' i; -\, 



(27T) 



(1) 



where 5(k) is the Fourier transform of the density pertur- 
bation 5(2) — [p(2) — Pm]/pm, p(2) is the local density at 
comoving position 2, p m is the background matter density 
and W(k, R) is a filter function in Fourier space. It is com- 
mon to relate the smoothing scale R to the corresponding 
variance of the linear density field 



o\R) = S(R)= l^ k ^iw(k,R).. 



k 2tt 2 



(2) 



where P(k) is the matter power spectrum in linear perturba- 
tion theory. We can refer to a trajectory S(x, S) as a sequence 
of overdensities given by subsequent increases in the smooth- 
ing scale by increments AS. When a tophat filter in fc-space 
is used then 8(2, S) executes a random walk. Given an un- 
derlying Gaussian distribution for the linear density field, 
the excursion set formalism allows us to associate proba- 
bilities to random walks that satisfy a given set of criteria 
for the smoothing scale at which they cross various density 
thresholds. Its use in defining the statistics of objects in the 
nonlinear regime requires a model that associates such cri- 
teria to objects. 



2.2 Spherical evolution and SVdW model 

The spherical evolution model provides a complete descrip- 
tion of the nonlinear evolution of a spherically symmetric 
top-hat density perturbation. One of the main features of 
this model is that the evolution does not depend on the ini- 
tial size of the region, i.e. on the initial radius or enclosed 
mass, but only on the amplitude of the initial top-hat over- 
density. 

For the collapse of perturbations, the spherical evolu- 
tion model in combination with the excursion set provides a 
good description of the statistics of dark matter halos. As we 
review in Appendix [XJ collapse occurs when the linear den- 
sity fluctuation reaches a critical value or barrier S c . We can 



then use the excursion set formalism to determine the frac- 
tion of trajectories that cross this barrier for the first time, 
accounting for the cloud- in-cloud process, within some diner 
of a smoothing scale a through the differential fraction 



/in a (cr) 



d/ 
dln<7 



2 5 c 



(•3) 



Since both mass and number are conserved in the collapse, 
the linear theory mapping a(M) carries over to the nonlinear 
regime and so the mass function, or the comoving differential 
number density of halos is 

dn _ p m f ^ dmo-- 1 
dlnAf ~ M /Wl ' dlnM ' () 

where Pm/M is the number density of such objects if the 
fraction were unity. 

We can extend the model to underdense regions in the 
initial density field. These are naturally associated with 
voids in the evolved density field today. A key assumption 
in making the connection between the excursion set and the 
abundance of nonlinear objects is that each collapse occurs 
in isolation. This makes sense for collapsing objects since the 
comoving volume occupied shrinks. In contrast to overdense 
regions which contract, voids expand. We shall see that this 
causes a problem for mapping excursion set predictions onto 
the statistics of voids. 

Nonetheless let us s tart with the simple spher ic al evo - 
lution model following ISheth fc van de Weygaertl l|2004) . 
The critical density threshold is defined to be when the 
expanding shells cross (see e.g. Suto. Sato fc Satol 11984 ; 



iFillmore fc Goldreichlll984l : lBertschingerlll985h . As shown in 
Appendix |A1 for an Einstein de-Sitter (EdS) universe, this 
occurs when the nonlinear average density within the void 
reaches p v = 0.2p m or when the linear density threshold 
reaches S v — —2.7. Note we will use this notation of p v to 
refer to the nonlinear density of the void region and 8 V to 
refer to the linear underdensity used as a threshold in the ex- 
cursion set model. We show in Appendix[A] that these EdS 
values suffice for the accuracy to which we wish to describe 
alternate cosmologies such as the ACDM model. 

Once we have this value for the void barrier we can fol- 
low the excursion set formalism for determining the fraction 
of random walks which pierce the barrier S v . Similar to the 
cloud-in-cloud process, the void-in-void process accounts for 
the fact that a void of a given size may be embedded in an- 
other underdense region on a larger scale. We thus define the 
first crossing distribution by associating the random walks 
with the smoothing scale for which they first cross the bar- 
rier 5 V . 

The second process, the void-in-cloud scenario, occurs 
when a void of a given size is embedded in an overdense re- 
gion on a larger scale, which will eventually collapse to a halo 
and squash the void out o f existence. In order to account for 
the void-in-cloud effect, ISheth fc van de WeygaertJ (J2004h 
proposed that the excursion set method applied to voids 
requires a second barrier, the threshold for collapse of over- 
dens e regions, S c . In calculating the fi rst crossing distribu- 
tion, ISheth fc van de Weygaertl (J200J) argued that we need 
to determine the largest scale at which a trajectory crosses 
the barrier S v given that it has not crossed 5 C on any larger 
scale. They posit that the value of 5 C should lie somewhere 
in between S c — 1.06, the value at turnaround in the spheri- 
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Vdn 

linear k\\\\\\\N 
SVdW 




r [Mpc/h] 

Figure 1. Void abundance model predictions. In the SVdW 
model, the number density of linear underdensities (blue curve) 
remains unchanged in void formation and only their sizes change 
(arrow to orange curve). In the Vdn model, the number density 
also changes so as to conserve the volume fraction in voids, low- 
ering the amplitude at fixed shape (arrow to grey curve) . Varying 
1.06 ^ 8 C ^ 1.686 (shaded or hatched regions) changes the abun- 
dance significantly only for small voids r < lh~ 1 Mpc. We take 
5 V = —2.7 throughout. We use the erg = 0.8 ACDM cosmology as 
listed in Tablc[2]here and in the following figures unless otherwise 
stated. 



cal collapse mod el, and S c = 1.686, the value at the point of 
collapse (see also lParaniape. Lam fc Shethl2012 | ). In this pa- 
er we shall refer to the model of lSheth fc van de WevgaertJ 



J2004 ) as the 'SVdW model 

By the same reasoning as applied to halos, the SVdW 
formula for the abundance of voids in linear theory is given 
by 



dn L 
dlnAf 
where 



p m j / ^ dlna 

M fl**W dlnM 



jw(<r) = 2^V 



2 jttx sm(jirV), 



with 



V 



|<5v 



V 



(5) 



(6) 



(<) 



5c + |<M' \8 V \ ' 

Note that ISheth fc van de WeygaertJ (|2004T ) give j\ aS = 
Sdf/dS = fi-na/2. We have added the subscript "L" to re- 
mind the reader that the logic relies on equating a number 
density derived from linear theory to the number density of 
some nonlinear object for reasons that will be clear below. 
Since the infinite series in equation (J6]) is cumbersome 
to work with, it is useful to have an accurate closed form ex- 
pression. As we discuss in Appendix iBl the accuracy of the 
approximation given in lSheth fc van de Weygaertl ((200J) is 
uncontrolled asa-)oo. Instead we find the limiting forms 
for equation J6| such that the domain of validity of the ap- 
proximation is well-defined. Note that as a — > 0, the oppos- 



ing barriers are high and the sum must return the single 
barrier expression since the probability of first crossing the 
collapse barrier is vanishingly small. This fixes the form as 
x — > 0. As T> increases toward unity, we lower the collapse 
barrier relative to the void barrier and the value of x at 
which this limit is approached decreases. Correspondingly 
to achieve a matching at this point, we need to keep more 
terms in the sum. The largest value that we will be inter- 
ested in is T> < 3/4 and so it suffices to keep 4 terms 



/lna(cr) 



2 [<Jv| 



x < 0.276 



2 ELi e~^~ jnx 2 sin(JTvV), x > 0.276 



(8) 



which is accurate at the 0.2% level or better across the do- 
main of validity. The approximation of equation © is used 
in all the numerical work throughout the paper. 

We can alternately express the number density in terms 
of the linear theory radius of the void r^. Using p m /M = 
1/V(j"l) and defining the volume of a spherical region of an 
arbitrary radius, R, as 

V{R) = ^R\ (9) 

we obtain 

dn L _ /in<r(o-) din a -1 
dlnr L ~ V(r L ) dlnr L ' ( ' 

In the spherical evolution model, the actual void ex- 
pands from its linear radius. At the epoch of shell crossing 
p v — 0.2p m . Given that 

1/3 

(11) 






spherical expansion predicts that this expansion factor is 
r ~ 1.7rL- The void abundance therefore becomes 



dn 



dn L 



dlnr dlnr-L \r L =r/i.7 



(12) 



Note that in this model dn/d In r shifts left to right in scale 
through the nonlinear growth but does not change in ampli- 
tude, as is shown in Fig. [T] 

The SVdW model has two parameters S c and <5 V . The 
latter is fixed by the shell-crossing criterion whereas the 
former is expected to vary within 1.06 ^ 5 C sj 1.686. In 
Fig. [T] we also show that for the range of radius of interest 
(r > lh~ 1 Mpc), changing 8 C within its expected range has 
little effect on the void abundance. 

The SVdW model makes a very specific prediction for 
the abundance of large voids. Again the key assumption of 
the SVdW model is that the comoving number density of 
objects is conserved during the evolution n — til and only 
their size has changed. Unfortunately, for spherical evolution 
this assumption is invalid for large voids. In particular, the 
cumulative volume fraction in voids larger than R defined 



T(R) 



dr dn 

— v ( r )-n — . 

r d In r 



(13) 



exceeds unity for radii of interest. In Fig. [2] we demon- 
strate that this problem cannot be cured by changing 5 C 
within the expected range as it only affects small voids 
whereas the problem appears at R « 2h~ Mpc. Indeed if 
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SVdW mm 
Vdn 
linear k\\\\\\\N 




R [Mpc/h] 

Figure 2. The cumulative volume fraction in voids with radii 
larger than R for the various models: linear theory (blue striped 
region, R = tl), SVdW model (orange striped region, R = r), 
Vdn model (grey shaded region, R = r). Regions correspond to 
the expected range of 1.06 ^ <5 C ^ 1.686 and we take 8 V = —2.7 
throughout. For SVdW the fraction unphysically exceeds unity 
at R ?s 2/i — 1 Mpc while for the Vdn model conserves the total 
fraction from the linear theory of -F(0) ~ 0.3. 



we take R — » then for the exact /i n a given by equation ([6]) 
l|Sheth fc van de Wevgaertll2004f ) 



m=i&) r^/w = (- 



''L 



v-l 



[1-V). 



(14) 



This result suggests that reducing 5 V — > simultane- 
ously takes r — > ri, and T> —¥ bringing SVdW asymptoti- 
cally to physicality J-(Q) — 1. Strictly speaking, S v is fixed by 
the shell-crossing criterion. However, given the approximate 
nature of the correspondence between the isolated spherical 
expansion model and real voids, it is interesting to explore 
whether modifications to this criterion can bring the SVdW 
model into agreement with physicality and simulations. If 
we change the nonlinear density at which voids are defined 
Pv/pm, the linear density threshold 5 V and the expansion fac- 
tor r/rL must change in a self-consistent fashion (see Fig. lAll 
and equation (JBTJ). In Fig. [3] we show that changing 8 V alters 
the shape of the abundance function. As \8 V \ decreases, the 
steepness of the abundance function also decreases. Thus, 
although lowering 8 V can make the total volume fraction 
physical (Fig. [3] lower panel), it increases the abundance 
of the largest voids. We shall see that correspondingly the 
SVdW model greatly overpredicts the abundance of large 
voids regardless of the choice of 5 C and 8 V . 



2.3 Volume conserving Vdn model 

We propose a simple fix to the unphysicality of the SVdW 
model. We require that the volume fraction and shape of 
the abundance function is fixed during the expansion, rather 
than assuming that the expansion of isolated voids preserves 
their total number density. Specifically, if we define the vol- 
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Figure 3. Relaxing the shell crossing criteria of the SVdW model 
void abundance predictions. Upper: Variation of the void undcr- 
density p v /p m changes both the shape of the abundance through 
the linear barrier <5 V and the size of the voids or horizontal shift 
through r/rx = (Pv/pm) -1 ' 3 - Decreasing |<5 V | increases the num- 
ber of large voids and decreases that of small voids. Lower: The 
cumulative volume fraction in voids with radii larger than R de- 
creases as |<5 V | — > and R — >• but at the expense of making the 
larger voids more abundant. We use 8 C = 1.06 throughout. 



ume fraction in linear theory, J-l, as 
J Rlj r L dlnr-L 



(15) 



then this fraction is conserved if we define the nonlinear 
abundance as 



V(r)dn = V> L )dn L | rL( 



r) 



(16) 



In this picture, when a void expands from r^ — > r it combines 
with its neighbours to conserve volume and not number. 
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(17) 



Thus the abundance becomes 

An V{tl) driL din 7X1 

dlnr V(r) dlnrL dlnr lr L (V) 

flnaja) dlno- -1 dlnr L I 
V(r) dlnrL dlnr lr L (r) 

We call this model the VAn model and show its abundance 
prediction in Fig.Q] We have left the mapping r(rh) general 
here since the specific form from isolated spherical expansion 
until shell crossing may not apply here. We will however 
adopt r = 1.7rL for voids with nonlinear density p v — 0.2p m 
from N-body simulations as a starting point. Note that in 
this case dlnrL /dlnr = 1 and the impact of going from the 
linear to the nonlinear abundance is both a shift in scale and 
a change in amplitude with no change in shape, as is shown 
as the combination of arrows in Fig. Q] 

In Fig. [2l we also show the cumulative volume fraction 
with this abundance function, along with that for linear the- 
ory defined in equation (|15|) . Since by construction the vol- 
ume fraction is conserved, the two curves differ only by a 
horizontal shift in scale. 

Since the VAn model is not the unique means of con- 
structing a physical model, it is interesting to explore other 
ways of keeping the volume fraction below unity. Phe- 
nomenologically, we can decouple the relationship in equa- 
tion (|B1|) between the parameters S v and r/r-L provided by 
the spherical expansion model. In fact, we can choose these 
parameters so as to mimic the VAn predictions for a fixed 
cosmology. For example, in the upper panel of Fig. [4j we find 
we can change the parameters 5 V — > —2 and r/r-L — > 1 in the 
SVdW model to fit the VAn model in the a 8 = 0.8 ACDM 
cosmology listed in Table [2] However, this change then pre- 
dicts very different abundances than the VAn model for a 
different cosmology as shown with the EdS cosmology listed 
in Table [2] and Fig. [4] (lower panel) . We shall show below 
that simulation results favor the VAn model over universal 
changes in S v and r/rx,. 



mod SVdW, 8 v =-2.0, r/r L =1 .0 
Vdn, 8 v =-2.7, r/r L =1 .7 




mod SVdW, 8 v =-2.0, r/r L =1 .0 
EdS Vdn, 8 v =-2.7, r/r L =1 .7 




r [Mpc/h] 

Figure 4. Void abundance in the Vdn model (black dot-dashed 
curve) and a modified SVdW model (blue solid line) with ad hoc 
variations designed to fit the ACDM Vdn model. Upper panel: we 
choose S v = —2 and r/r^ = 1 in the SVdW model in violation of 
spherical expansion predictions in a ag = 0.8 ACDM cosmology. 
Lower panel: we show that the same set of parameters give a 
poor fit, in an EdS model (sec Table [2} • For all curves we use 
<5 C = 1.686. 



3 SIMPLE VOID FINDING ALGORITHM 



In Section T3. II we outline the void finding algorithm used to 
identify voids in both the dark matter and halo populations 
in this work. In Section f3. 21 we present the details of the N- 
body simulations which were carried out as summarised in 
Tabled 



3.1 Void finder 

As we have already mentioned one of the main complications 
in studying the distribution of voids in the large scale struc- 
ture of the U niverse, is finding a robust definition of what 
a void is fsee lColberg et alj|2008i . for a comparison of void 
finders). In this work we wish to make a direct comparison to 
the predictions of the excursion set formalism which assumes 
that these underdense regions are non-overlapping spheres 
of a given underdensity corresponding to a region at the mo- 
ment of shell crossing. We shall retain the moment of shell 
crossing as the key feature which defines a nonlinear void in 
the matter distribution today although we also compare the 
measured abundance of voids with different underdensities 



to the predictions of the volume conserving model in Section 

ED 

We start w ith the publicly available code ZOBOV 
(Ncyrinck 2008) which uses Voronoi tessellation to estimate 
densities and find both voids and subvoids. The main ad- 
vantage of using tessellation methods is that it gives a local 
density estimate by dividing space into cells, where the cell 
around any given particle is the region of space closer to 
that particle than to any other. The Voronoi tessellation 
also gives a natural set of neighbours for each particle which 
ZOBOV uses to construct zones around density minima. 

The output from ZOBOV is useful for our purposes for 
two main reasons. Firstly, it outputs a linked list of zones 
in the dark matter distribution, which is also ordered by 
density contrast. It thus provides a tree structure which we 
can prune according to the definition of a void. Secondly, 
ZOBOV identifies the 'core' or least dense particle in a zone 
and returns its density as well as a measure of the probabil- 
ities that each collection of zones arises from Poisson fluc- 
tuations. Note that the list which ZOBOV returns contains 
zones of various densities and Poisson probabilities, some of 
which could be overdense or not statistically significant, so 
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Table 1. A toy example showing the first five columns from a 
ZOBOV output file. 



Void# (zones) FileVoid# CorePar 
1 (a, b, c) 2945 26 

2 (b, c) 5033 83 

3 (c) 1814 45 



CorcDcn ZoneVol 

1.08e-02 4.9e+03 

1.8e-02 7.8e+02 

2.0e-01 1.9e+02 



it is necessary to prune the output from ZOBOV in order 
to construct a void catalogue. 

In constructing our void finder the goal is to identify 
all spherical non-overlapping underdense regions of average 
density p v = 0.2p m in a dark matter simulation. We use 
the output from ZOBOV and find spherical regions centred 
around the core particle (lowest density particle) in a zone, 
which can encompasses any particles which are around the 
zone returned by ZOBOV and not necessarily part of the 
particular zone or collection of zones returned by ZOBOV. 

One of the outputs from ZOBOV is a text file which 
lists individual zones and joined zones which are added 
to the list in a sequential process analogous to water 
flooding a plane with troughs of variou s heights (see also 
IPlaten. van de Weygaert fc Jones! l2007f ). During flooding, 
when water from a particular zone or joined zones flows 
into a neighbouring deeper zone, the process stops and the 
zone is recorded in the list. A toy example of the output is 
shown in Table [T] where the zones are listed in order of den- 
sity contrast. Here FileVoid# and CorePar refer to unique 
identification tags for the void and its core particle respec- 
tively. CoreDen is the density, in units of the mean, of the 
void's core particle. 

In order to count non-overlapping regions with an av- 
erage density of 0.2, we perform the following two stages of 
analysis on the output from ZOBOV. Firstly, starting from 
the top of the ZOBOV output file we determine if the first 
collection of zones listed e.g Void #1, which is made up of 
zones a,b and c, pass the following criteria: 

• [r m i n ,r max ]: The radius corresponding to a sphere of 
equal volume should be > r m i n and < r max . 

• The core particle density is < 0.2p m . 

As we are searching for regions which have an average den- 
sity p v = 0.2p m we also only consider zones in the list which 
have a density p v ^ 0.2p m in order to speed up the search. 

If the collection of zones fulfils all of the above, then 
we proceed to the second stage. A spherical region around 
the centre is found by iteratively including one particle at 
a time moving away from the centre of the void until p v = 
0.2p m . We assume that the volume corresponds to a sphere 
with radius equal to the distance from the centre to the last 
particle included. If for example, Void #1 did not pass the 
r max , r m in criteria then we consider if the deepest zone, zone 
o, does and if so we grow a sphere around the centre of 
this single zone. This step is important as the output from 
ZOBOV only lists the deepest zone 'a' once and if Void 
#1 fails the r max ,r m i n criteria the void finder would miss 
counting the deepest zone in the simulation box. 

We then proceed to the next line in the output file and 
perform the same two stages of analysis. All of the particles 
in the spherical regions which are grown have been tagged 
and at any stage if there is any overlap of spheres we dis- 
regard the less underdense zone to avoid double counting 



any volume in the simulation. Note if we consider the out- 
put from ZOBOV as a tree structure then this procedure 
is similar to walking the tree from root to tip, pruning any 
branches after our criteria are met. 

Using a cut in r max and r m i n as above allows us to avoid 
considering spuriously small voids and the first output in the 
text file which is a void which takes up the entire simulation 
box, however we have checked that our results are not sensi- 
tive to changes in r max and r m i n but we retain these criteria 
in the void finder to speed up computation. For our simu- 
lations we use the following: [r m i n ,r max ] = [0.5, 15], [1, 30], 
[2,60] and [4,120] /i _1 Mpc for the 64/i _1 Mpc, 128/i _1 Mpc 
and 256/i _1 Mpc and 500/i _1 Mpc boxes respectively. 

We tested several different criteria in identifying voids 
and found that the two points listed above are sufficient to 
identify significant non-overlapping regions of a given under- 
density in the particle distribution. In testing the robustness 
of our void finder we considered the impact of the following 
adjustments to the method: 

• As an alternative to using the core particle to define 
the centre of the void we can use a volume-weighted centre 
defined as 






(18) 



where xl and Vi are the position and volume of each parti- 
cle in the zones returned by ZOBOV respectively. We found 
that the abundance of voids was not substantially affected 
by this choice and so we use the core particle as the centre 
of the spherical region. Note that using the volume- weighted 
centre is more robust w hen using stacked voids (see e.g. 
lLavaux fc Wandelt|[2012T ). 

• Instead of allowing the spherical region to include all 
particles which are around the void we consider restrict- 
ing it so that only particles which ZOBOV list as being 
part of the zone are included when growing the sphere. We 
find that in the majority of cases the region of underden- 
sity p v = 0.2p m is contained within the collection of zones 
returned by ZOBOV and so this restriction does not affect 
the measured abundance of voids. Our results in Section [3] 
use all particles within 1.5 times the radius of the zones to 
find the spherical void. Including particles only within this 
radius was found to be sufficient considering the original 
collection of zones was required to have an average den- 
sity of > 0.2p m . An alternative approach to this would be 
to use the actual volume of each zone particle when trying 
to find a void of a given average density. This would allow 
for irregularly shaped voids which it could be argued is a 
more 'natural' description of an actual void, however as we 
mentioned we are trying to compare with the excursion set 
theory for abundances which assumes spherical voids. 

• We originally included a third criteria in our void finder 
by requiring that the probability a zone arose from a Poisson 
process was less than a given significance (see lNevrinclJ2008l . 
for more details). However we found that the core particle 
density requirement by itself was sufficient to get r id of spu- 
rious voids. This also agrees with the findings of iNevrinckl 
(|2008h . 

• In the algorithm we have described we stop growing a 
sphere around the core particle when we find the desired 
underdensity at the maximum radius at which this occurs 
within the radius of the collection of zones. This is a different 
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Table 2. Details of the simulations used in this work. 



Model 


i l m 


h 


CT8 


Boxsize 


# particles 


# 


simulations 


z output 


ACDM 


0.26 


0.715 


0.8 


64, 128, 256 
500 
500 


256 3 
400 3 
256 3 




8 
1 
8 


0, 0.5, 1 





ACDM 


0.26 


0.715 


0.9 


64, 128, 256 


256 3 




8 


0, 0.5, 1 



EdS 



0.7 



0.8 64, 128, 256 



256 3 



approach to simply stopping to record the first radius where 
p v = 0.2p m which would not take into account void-in-void 
scenarios. In practice we find that accounting for a void- 
in-void effect alters the measured abundances by a small 
amount (e.g ~ 7% over the range 1 < r(/i Mpc) < 10). 

• The above method does not allow any overlap of voids 
within the simulation in order to compare with the excursion 
set method. In practise for regions of p v = 0.2p m we found 
the overlap was very small for the simulations we consider. 

• ZOBOV is run using all the particles in the simulation 
with a run-time density threshold parameter which can limit 
the growth of a collection of zones into high-density regions. 
We set this parameter to 0.2 however we have verified that 
changing this run time parameter has little effect on the 
abundance of underdensities found by our void finder. 

3.2 TV-body simulations 

We measure the abundance of voids in the dark matter 
distribution using a series of N-body simulations in vari- 
ous box sizes. These simulations were carried out at the 
University of Chicago u sing the TreePM simulation code 
Gadget-2 (jSpringefecioBT ). The ACDM model used has the 
following cosmological parameters: fi m = 0.26, Qde = 0.74, 
Q h = 0.044, h = .715 and a spectral tilt of n s = 0.96 
(jSanchez et al.ll2009r i. The linear theory rms fluctuation in 
spheres of radius 8 ft -1 Mpc is set to be erg = 0.8 for our main 
simulation set of 8 independent realisations of the ACDM 
cosmology. In order to investigate the abundance of voids in 
different cosmologies we also carry out two additional sim- 
ulations; one with a ACDM cosmology and erg = 0.9 and 
another which we refer to as the 'EdS' simulation which has 
fi m = 1. The EdS simulation is not a viable cosmological 
model for our Universe as it has already been ruled out by 
many observations but we use it here as a tool to exam- 
ine how robust our void models are to large changes in the 
power spectrum or cosmology. 

The simulation details are summarised in Table [2] Most 
of the simulations use TV = 256 3 particles to represent 
the dark matter while for the larger simulation box of 
500ft _1 Mpc we use 400 3 particles. The error on the abun- 
dance of voids measured in the 500ft" 1 Mpc box is estimated 
from eight lower resolution simulations which have 256 par- 
ticles in a computational box of 500/i _1 Mpc on a side. These 
lower resolution simulations have a mean abundance which 
agrees with the 400 particle simulation over the range of 
scales which we consider and are computational less expen- 
sive to run and analyse with the void finder. The initial con- 
ditions of the partic le load were set up with a glass config - 
uration of particles (jBaugh. Gaztanaga fc Efstathioulll995r ) 
and the Zeldovich approximation to displace the particles 
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Figure 5. Void abundance in simulations vs. predictions with 
p v = 0.2p m in the dark matter distribution of the erg = 0.8 ACDM 
cosmology in simulation box sizes 64h -1 Mpc (green), 128fe _1 Mpc 
(purple), 256/i — x Mpc (red) and 500/i _1 Mpc (cyan) on a side. The 
error bars represent the scatter on the mean from eight differ- 
ent realisations of this cosmology in each box size. The range in 
predictions cover the parameter interval 5 C = [1.06,1.686] with 
<5 V = —2.7 and are consistent with simulations for Vdn (grey 
shaded) but not SVdW models (orange hatched). 



from their initial positions. We chose a starting redshift of 
z = 100 in order to li mit the discreteness effects of the initial 
displacement scheme (jSmith et al.ll2003l ). The linear theory 
power spectrum used to generate th e initial conditions wa s 
created using the CAMB package of iLewis fc Bridle! (|2002h . 
Snapshot outputs of the dark matter distribution as well as 
the group catalogues were made at redshifts 1, 0.5 and 0. 
In the following section we also test voids that are identi- 
fied with dark matter halos. The simulation code Gadget2 
has an inbuilt friends-of-friends (FOF) halo finder which was 
applied to produce halo catalogues of dark matter particles 
with 10 or more particles. A linking length of 0.2 times the 
mean interparticle separation was used in the halo finder. 



4 RESULTS 

In the following sections we compare simulation results for 
the abundance of voids with the predictions of both the vol- 
ume conserving and the SVdW model. In Section 14.11 we 
present the measured abundances of voids from the ACDM 
simulations in four different simulation box sizes. In Sections 
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Figure 6. Void abundance for different defining underdensities p v = 0.3p m (left panel), p v = 0.4p m (right panel). The grey shaded 
region represents the excursion set predictions with varying amplitude and using a linear underdensity value <5 C = 1.686 and <5 V , given in 
the legend in each panel. The amplitude rescaling vs the SVdW predictions ranges from p v /pm (Vdn; top black dashed curve) to 1/5 
(bottom black dotted curve) which both preserve agreement for p v = 0.2p m - 



14. 2114. 3114. 41 we test the robustness of the models to variation 
in the critical void underdensity, redshift, and cosmology re- 
spectively. In Section 23] we present the abundance of voids 
identified in the dark matter halo catalogue. 



4.1 Baseline model comparison 

We implement the void finder, which is described in Sec- 
tion [3TT] to measure the abundance of spherical voids which 
have p v = 0.2p m at z = in all four simulation box sizes of 
the ACDM cosmology, see Table[2] Fig.[5]shows the average 
number density as a function of radius, of voids measured 
from eight different realisations of the ACDM cosmology in 
simulation box sizes 64/i~ 1 Mpc (green), 128/i _1 Mpc (pur- 
ple), 256/i _1 Mpc (red) and 500/i -1 Mpc (cyan) on a side. 
The error bars represent the scatter amongst these simula- 
tions. 

The orange hatched region in this figure represents the 
SVdW model within the parameter interval S c — [1.06,1.686] 
and S v = —2.7 and assuming that the voids have expanded 
by a factor of 1.7 today. The grey shaded region shows the 
Vdn model for the same parameters. As discussed in Section 
12.21 the range in S c accounts for the void-in-cloud process by 
which a void in a larger overdense regions will be crushed out 
of existence. As we can see from Fig. [S] this only affects the 
smallest voids of r < l/l - Mpc and for larger voids the abun- 
dance is insensitive to 8 C . The decrease in the void abun- 
dance at r(/i _1 Mpc) ~ 2.5, 1.5 and 1 for the 256h _1 Mpc, 
128ft _1 Mpc and 64ft -1 Mpc boxes shows the resolution limit 
for each of these simulations where small voids are not fully 
resolved and so the abundance is decreased. 

The SVdW model overpredicts the abundance by a fac- 
tor of 5 whereas the Vdn model agrees with simulations 
to ~ 16% across the range 1 < r(/i _1 Mpc) < 15 where 
the results measured from simulations in different box sizes 
has converged. The Vdn model conserves the volume rather 
than the number of voids and hence implies that the number 



density decreases in going from the linear to the nonlinear 
regime by the same amount that the volume of the voids 
grow. It is somewhat surprising that using the factor of 1.7 
in this model, which applies to the expansion of isolated ob- 
jects, fits the results from the simulations where voids have 
merged as they expand. It is important to test that this is 
not just a coincidence but rather is robust to other choices 
of parameters in the simulations. 



4.2 Underdensity variation 

In both the SVdW and Vdn models, we adopt the shell 
crossing criteria p v = 0.2p m for defining the void and match 
predictions to p v as denned by the simple void finder of 
Section 13.11 If the agreement between the Vdn model and 
simulations was robust, we would expect that it would be 
preserved for at least small variations in this criteria. 

We modify our void finder such that the largest non- 
overlapping spherical regions which have densities p v = 
0.3p m and p v = 0.4p m are recovered from the simulations. 
The results are shown in the left and right panels of Fig. [5] 
respectively. The errors plotted in this figure represent the 
scatter on the mean from eight simulations. 

As discussed in Section \2. 21 (see also Fig. lAll and equa- 
tion (|B1|| ). changing the underdensity criteria in the spher- 
ical evolution model alters the shape of the abundance 
function through the linear threshold 8 V . Specifically, for 
p v = 0.3p m , S v — —1.8; while for p v = 0.4p m , 5 V — —1.24. 

For dark matter voids with p v = 0.2p m the predicted 
abundance for the Vdn model are approximately a factor 
of 5 smaller then those of the SVdW model. In modelling 
the number density of underdense regions with p v = 0.3p m 
and p v = 0.4p m , which cannot be directly related to shell 
crossing in the spherical expansion model, we adopt a phe- 
nomenological approach. Given that for p v = 0.2p m , the 
Vdn model has the same shape as the SVdW model but a 
factor of Pv/pm = 1/5 lower amplitude, we can preserve the 
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good match there by either following the Vdn prescription 
literally and rescaling SVdW by p v /p m or by simply keeping 
this factor fixed at 1/5. 

This range is plotted in Fig. [5] as grey shaded regions 
bounded by the two limiting cases, black dashed and dotted 
lines. Simulation results clearly favor the simple phenomeno- 
logical prescription of rescaling the amplitude by 1/5. The 
Pv/pm scaling prescription of the Vdn model would over- 
predict the amplitude by approximately 1.5 for voids with 
p v = 0.3p m and 2 for p v = 0.4p m . These results again high- 
light the point that excursion set models predict the overall 
shape of abundance function accurately and only the am- 
plitude needs to be altered to fit the simulations results, 
here without the benefit of volume conservation as motiva- 
tion. Note that the preferred rescaling of 1/5 is more than 
sufficient to bring the predictions to a physical void filling 
fraction for p v ^ 0.2p m . 



4.3 Redshift variation 

Next we check the robustness of results to the redshift at 
which the void abundance is measured. Fig. shows the 
number density of voids as a function of radius at z = 0.5 
(blue) and z = 1 (red) measured from the ACDM, as = 0.8 
simulation. The measured abundances from the three simu- 
lation box sizes 64/i _1 Mpc, 128/i~ 1 Mpc and 256/i~ 1 Mpc are 
the volume-weighted averages and errors over 8 realisations. 
The volume conserving (Vdn) model is shown as a black 
hatched (grey shaded) region using <5 V = —2.7 at z — 0.5 
(z = 1) and the parameter range <5 C = [1.06, 1.686]. Note 
we have used only one colour for the results from the three 
simulation boxes at each redshift for clarity in this figure. 

We again find that the Vdn model works very well in 
reproducing the abundance of voids in the dark matter in a 
ACDM universe at both redshifts, while the SVdW model, 
which is not plotted here for clarity, again overpredicts the 
abundances by approximately a factor of 5. Fig. [7] shows 
that smaller (larger) voids are more (less) abundant at z — 
1 compared to z = 0.5 which is also found in the model 
predictions at both redshifts. 



4.4 Cosmological parameter variation 

In order to check if the Vdn model for the abundance of 
voids works when we change the cosmological model we 
have run two simulations of alternative cosmologies to the 
standard ACDM with erg = 0.8 which we discussed in the 
previous section. In the first alternative cosmology we have 
chosen to modify only the value of erg to 0.9, see Table[2] our 
second alternative cosmology is an Einstein de-Sitter (EdS) 
universe where the matter density parameter Jl m = 1. The 
linear perturbation theory power sp ectra for these simula - 
tions were generated using CAMB (jLewis fc Bridle! 120021 ) 
and normalised to erg = 0.9 (erg = 0.8) for the ACDM (EdS) 
simulations in order to generate the initial conditions for 
the simulations and the variance a(R) which is used in the 
excursion set model for the abundance. 

The measured z — number density of voids with p v = 
0.2p m in the erg = 0.9 and EdS simulations are shown in Fig. 
[8] The volume conserving model is shown in both panels as 
a grey shaded region as in previous plots. We have used the 
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Figure 7. Redshift dependence of the void abundance with 
p v = 0.2pm at z = 0.5 (blue) and 2 = 1 (red) measured 
from the ACDM, erg = 0.8 simulation. The black hatched (grey 
shaded) region represent the Vdn model using linear underden- 
sity values of <5 V = —2.7 at z = 0.5 (z = 1) for the range 
(5c = [1.06,1.686]. Note the measured average abundances and 
errors from the 256h _1 Mpc, 128/i _1 Mpc and 64h _1 Mpc simu- 
lation boxes are the volume-weighted values. Note in this figure 
we have plotted the results from the three simulation boxes using 
the same colour for clarity. 



same value, 8 V = —2.7, for the linear perturbation theory 
underdensity. Note this parameter is different in different 
cosmologies, however we find that such a small change in 5 V 
going from an EdS to a ACDM universe has a small impact 
on the predicted abundance of voids in the excursion set 
theory and the main differences arise from the change in the 
variance, a(R) (see Appendix |X) |. 

From Fig. [8] it is clear that the volume conserving 
model works well in both of these cosmologies and fit the 
abundance of voids to within 25% over the range 1 < 
r(h~ Mpc) < 15. It is interesting to note the overall de- 
crease in the abundance of voids in the dark matter distri- 
bution for voids with small radii r < 2h~ Mpc in these two 
cosmologies which is most obvious in the measured number 
density from the EdS simulation and a larger abundance for 
the erg — 0.9 cosmology for large r. It is also clear from Fig. 
[8] (lower) that the excursion set model predicts more squash- 
ing of smaller voids due to the void-in-cloud effect but this 
is occurring right on the resolution limit of our simulations 
at 2 < r(/i _1 Mpc). Finally note that even if we modified the 
SVdW model in the ad hoc manner of Fig. [4] to match the 
simulation results of ACDM with erg = 0.8, the predictions 
would be far off simulation results for the EdS cosmology. 



4.5 Halo denned voids 

Voids in the galaxy population are defined not through the 
dark matter density field but by the number density field rih 
of the dark matter halos they populate. In this section we use 
density minima in the halo number density. Our goal is to 
test how faithfully the abundance of voids in the dark matter 
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Figure 8. Void abundance for alternate cosmological parameters 
at z = 0. Upper: ACDM with the initial conditions normalised to 
give erg = 0.9. Lower: EdS model with erg = 0.8 with f2 m = 1. The 
grey shaded region shows the Vdn model within the parameter 
interval <5 C = [1.06,1.686] and using S v = —2.7. 



matches that in the halo populations within the context of 
the simple void finder of Section [3TT] 

In this section we use the FOF halo catalogues from 
the 128, 256, 500 h~ Mpc simulation boxes and the pub- 
licly available hal o catalogues from the MultiDark and Bol- 
shoi simulations (|Riebe et al.l |201lj) which have computa- 
tional box sizes of L - 1000/i _1 Mpc and L = 250/i _1 Mpc 
on a side respectively. These halos have been identified using 
the B ound-Density-Maxima algorithm (JKlypin fc Holtzmanl 
1 19971 ). We only use halos which have Knax > 200km/s and 
M > 1 - 2 x 10 12 /i _1 Af Q from the Bolshoi and MultiDark 
simulations to ensure that the statistics are robust. We use 
the void finder described in Section [3. II to identify voids in 
the distribution of halos which have n v = 0.2nh where n v 
is the average number density in the void whereas n^ is the 
average in the whole simulation. Our final sample consists of 
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Figure 9. The number density of voids with n v = 0.2n^ in the 
halo distribution from the 128 (purple), 256 (red), 500 (cyan) 
h~ 1 Mpc simulation boxes. The results from the MultiDark (Bol- 
shoi) simulation are shown in dark green (blue). The error bars 
represent the error on the mean from eight simulations. The er- 
rors on the MultiDark simulation represent the Jackknifc error on 
the mean. The grey shaded region bounded by the black dashed 
and dotted line represents the volume conserving model with 
<5 V = —1.24 and varying amplitude as in Fig. [6] The grey solid 
line represents the Vdn model with <5 V = —2.7. 



5,768 voids using 1.7 xlO 6 halos from the MultiDark simu- 
lation and 4,826 voids using 2.2 x 10 6 halos from the Bolshoi 
simulation. Both of these simulations are of a higher resolu- 
tion than the ones we carried out in 128, 256, 500 /i -1 Mpc 
simulation boxes - it is useful to compare the abundance 
of voids in the halo population from these simulations to 
ours as an indication of the scales at which our results have 
converged. 

The measured abundance of voids in the halos popula- 
tion from our simulations and the Bolshoi and MultiDark 
catalogues at z = are shown in Fig. [9] The errors shown 
on the results from the 128 (purple), 256 (red), 500 (cyan) 
/i -1 Mpc simulation boxes are measured from the scatter 
amongst eight different realisations in each box size. The 
errors on the MultiDark simulation were obtained by jack- 
knife sampling from each simulation by dividing the sim- 
ulation volume into iV S ub = 8 equal subvolumes and then 
systematically omitting one subvolume at a time in order to 
calculate the void abundance on the remaining iV su b — 1 vol- 
ume. We find that voids identified in this manner through 
the halo distribution do not follow the Vdn model assum- 
ing n v = 0.2nh which corresponds to dark matter voids of 
p v = 0.2/5 m . They also do not follow the SVdW model which 
would have the same shape but 5 times the amplitude. 

Overall Fig.[9]shows that our void finder finds large halo 
defined voids that do not correspond to dark matter defined 
voids of the same underdensity for r > 10/i _1 Mpc. Although 
it is difficult to compare these results with previous work 
due to the large differences in the voi d finders used, quali- 
tatively this agrees with the findings of lBenson et al.l (|2003l ) 
who measured the void probability function from simula- 
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Figure 10. Left: A 10 X 50 X 50h~ x Mpc slice through the 500 h — 1 Mpc simulation box centred on a large, r ~ 21h~ 1 Mpc, void in the 
halo population (black dots). The diameter of the void is shown as a dashed grey line and the coloured contours represent the log of 
the densify field which has been evaluated on a grid of 256 3 points. The red circles represent all the voids in dark matter which have 
p v = 0.2p m and whose centres are within 10h~ 1 Mpc of the centre of the void in the halo distribution. Right: A 60 X 14 X 60/i _1 Mpc slice 
through the 500 h _1 Mpc simulation box centred on a large, r ~ 26h~ 1 Mpc, void in the halo population (black dots). The red circles 
represent all the voids in dark matter which have p v = 0.2p m and whose centres are within 10h _1 Mpc of the centre of the void in the 
halo distribution. Note these voids only appear to be overlapping due to the projection effect. 



tions and found that the VPF for voids with r > 5ft _1 Mpc 
was much higher for the galaxy catalogues compared to the 
dark matter. These results illustrate the fact that there is 
not always a 1:1 correspondence between voids in the dark 
matter and the dark matter halo distributions and this is es- 
pecially pronounced when we define a void as having a fixed 
underdensity which is the same for dark matter and halos. 

To illustrate the mismatch between the voids which we 
find in the dark matter and halo distributions using the same 
underdensity criterion, Fig. 1 101 shows a 10 x 50 x 50ft _1 Mpc 
(left panel) and a 60 X 14 X 60/i _1 Mpc (right panel) slice 
through the dark matter density field which has been eval- 
uated on a grid of 256 3 points from the 500 /i _1 Mpc sim- 
ulation box. The coloured contours represent the log of the 
density field in each cell and the halos around each void 
are represented by black dots. The radius of each void is 
r ~ 21ft~ 1 Mpc (left panel) and r ~ 26/i _1 Mpc (right panel) 
and is shown as a grey dashed line in Fig. 1101 The red cir- 
cles in this plot show the voids identified in the dark matter 
whose centres are within 10/i -1 Mpc of the centre of the void 
in the halo population. Not only is it possible to find more 
than one dark matter void which overlaps with the halo void 
but the radii of the dark matter voids at which p v — 0.2p m 
are a lot smaller than the halo voids which satisfy the anal- 
ogous criterion. 

There are at least two possible ways to reconcile the 
Vdn predictions for the abundance of dark matter voids with 
that of the halo voids. Firstly, a scale dependent modifica- 
tion to the barrier in the Vdn model could be used to alter 
the underdensity threshold used to find voids in the dark 
matter. Secondly if we keep a fixed underdensity threshold 
to define dark matter voids, it may be possible to find a 
scaling of this threshold to define voids in the halo distribu- 



tion. These ideas are beyo nd the scope of this work but see 
iFurlanetto fc Piranl (|2006l ) for related ideas. 

As a simpler illustration of these ideas, in Fig. [9] we 
also plot the Vdn model assuming that halo defined voids 
of n v = 0.2nh correspond to dark matter defined voids of 
p v = 0.4p m . These predictions are plotted as a grey shaded 
region allowing the amplitude to vary from the predictions 
of the Vdn model which rescales the SVdW amplitude by 
Pv/pm = 0.4 (black dashed lines) and the rescaling of 1/5 
that fits our dark matter voids well as in Fig. [S] (black dotted 
line). Compared to the predictions of the Vdn model for 
dark matter voids of n v = 0.2rih (solid grey line), these black 
dashed and dotted curves match the abundance of voids 
in the halo populations better though no single rescaling 
matches perfectly across the full range. 



5 SUMMARY AND CONCLUSIONS 

The ne xt generation of galax y redshift surveys such as Big - 
BOSS (|SchIegel et all 120091), Euclid jLaureiis et all I20I1I ) 
and WFIRST (|Albrecht et al.ll2006l : iGreen et al.ll2012l ) will 
allow us to study the large scale structure of our Universe 
in ever greater detail. Cosmic voids represent one of the 
main components which strongly influence the growth of 
clusters, walls and filaments in the mass distribution. Study- 
ing the statistics and dynamics of these underdense regions 
is a promising way to test the cosmological model and hier- 
archical structure formation. 

Several challenges which may affect the usefulness of 
voids as a probe of cosmology are addressed in this pa- 
per, such as, given a survey or numerical simulation of a 
given size, how robust are the statistics on the number den- 
sity of voids of a given size, how accurately can we pre- 
dict the number density of these voids and how faithfully 
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do voids in the halo population represent voids in the dark 
matter. Using N-body simulations of a ACDM cosmology 
we test the excursion set model for the abundance of voids 
includ ing the model provided bv lSheth fc van de Weygaerj 
1)20041 ) . which takes into account the void- in- void and void- 
in-cloud scenari os. Our void finder makes use of the ZOBOV 
(Ncyrinck 2008) algorithm which uses Voronoi tessellation 
to locate density minima. We define a void as a spherical 
region around these minima with p v — 0.2p m and make 
use of several different computational box sizes so we can 
determine the volume and resolution that is needed in or- 
der to recover robust statistics for voids of a given size. We 
have tested the robustness of our void finder to the following 
changes and have found convergent results: using different 
simulation boxsizes and particle numbers, using a volume- 
weighted centre or the core particle to define the centre of a 
void, using all particles around the density minima or only 
particles in a zone given by ZOBOV and using only statis- 
tically significant voids or voids with a core particle density 
< 0.2p m . 

We find that the measured abundance of voids at z = 
fro m a ACDM simulation does not m atch the predictions of 
the ISheth fc van de Weygaerj (J2004J) model, which greatly 
overpredicts the results using a linear underdensity of 5 V — 
—2.7. Instead we find a volume conserving model, which is 
also based on the excursion set method with S v = —2.7, 
matches the measured abundances to within 16% for void 
radii 1 < r(h~ 1 Mpc) < 15. This model works remarkably 
well and suggests that the number density of voids decreases 
in going from the linear to the nonlinear regime by the 
same amount that the voids expand. This agreement is ro- 
bust to varying the redshift in the ACDM model as well 
as the underlying cosmology. Using simulations of different 
cosmological models, a ACDM cosmology with as — 0.9 
and a Einstein-deSitter cosmology with fi m = 1, we find 
that the volume conserving model works well and repro- 
duces the measured number density from each simulation 
to within 25% over the range 1 < r(/i _1 Mpc) < 15. We 
have also tested model predictions for density thresholds 
of p v = 0.3p m and p v = 0.4p m and find that the volume 
fraction physicality rescaling factor remains fixed at ~ 1/5 
rather than scaling as p v /p m - 

Using the number density threshold criteria of n v = 
0.2nh in our void finder we have examined the voids in the 
halo population from the 128, 256 and the 500/i -1 Mpc on 
a side computa tional boxes . We also use the Bolshoi and 
the Multidark (jRiebe et al.l l201lh simulations to measure 
void abundances. These two simulations are of a higher res- 
olution than our simulations and we have confirmed that 
our measured void abundances in the halo population agree 
with both the MultiDark and Bolshoi measurements. We 
find that the number density of voids in the dark matter 
halo distribution is very different to that in the dark mat- 
ter, with fewer small, r < 10/i _1 Mpc, voids and many more 
large, r > 10/i" 1 Mpc, voids. These results indicate that a 
given void in the halo distribution of fixed underdensity of 
n v = 0.2rih cannot be unambiguously related to a void in the 
dark matter of equal underdensity and in the case that there 
is a 1:1 correspondence, the radii at which a dark matter or 
halo void have a given underdensity can be very different. 

Cosmic voids are a promising and interesting tool which 
can be used to test many aspects of the ACDM cosmological 



model and having an accurate model for the number density 
of voids in the Universe represents a first step. In this work 
we have presented a model based on the excursion set the- 
ory which conserves the volume fraction of voids and works 
well at predicting the abundance of voids identified in dark 
matter from N-body simulations over a range of scales. Es- 
tablishing the relationship between voids in the dark matter 
and in the halo or galaxy distribution requires further study. 
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APPENDIX A: 
MODEL 



SPHERICAL EVOLUTION 



In this Appendix we review the spherical evolution 
model which describes the nonlinear evolution of an 
un-compensated spherically symmetric tophat underdense 
(overdense) perturbation. To illustrate the physics, we will 
use the spherical tophat model in Einstein-de Sitter cosmol- 
ogy as an example, which is analytically solvable until shell 
crossing. 

Consider an initial spherical tophat density perturba- 
tion (| Jo | <S 1) of physical radius Ro at ai = a(ti). We can 
think of the perturbation as composed of concentric mass 
shells, labeled by their initial physical radii Ri at ai. Let 
A(Ri,a) denote the average overdensity of the region en- 
closed by the mass shell Ri. Then its initial value is 



Ai(i?i) = A(Ri,a,i) 



So 
S (R /Ri) 3 



Ri < -Ro, 
ft > Ro- 



(Al) 



For brevity we omit the _Ri argument of Ai below. According 
to Birkhoff's theorem, the evolution of the mass shell _Ri 
only depends on the total mass inside Ri, but not the mass 
distribution or the mass outside. Thus the shell Ri evolves 
in the same way as a FLRW universe 



R(t;Ri 



R(t;Ri) 



H? 



(1 + Ai) 



-A, 



(A2) 



where -R(i) is the physical radius, and initial conditions are 
set to the growing mode in linear theory. By introducing the 
dimensionless conformal time 



Ri 



-A; Hidt, 



(A3) 



we can solve equation ()A2|) in a parametric form (to leading 
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order in Ai) 



R_ 
Ri 

Hit 



h 



f* 



j (cosh 77 — 1) 
1(1 — cos 77) 

"2 J (sinh 77 — 77) 
\ (77 -sin 77) 



So < 0, 
So >0; 

<5o < 0, 
So > 0. 



(A4) 
(A5) 



Given this evolution of mass shells, we are particularly 
interested in shell crossing for the expansion case (So < 0), 
which is usually seen as a characteristic event that signifies 
the formation of the void at a nonlinear level. Note that these 
solutions represent a family of trajectories labeled by Ri 
and parametrized by 77^ . We can find out when and where 
shell crossing first occurs by differentiating the parametrized 
solutions with respect to -Ri and 77, and requiring that dR 
and d£ vanish, for _Ri > Rq 



An A 12 
A21 A22 

where 

An = 2 



dRi/Ri 
d-q 



0. 



(A6) 



§* 



(cosh 77 — 1), 



1 5 a 

3 
9 5 5 

A21 = - -Ai (sinh?/ -77), 

_ 3 
15 5 

A22 = - -Ai (cosb.77-1). 



(A7) 



For this homogeneous system of linear equations to have 
nonzero solutions, we must have det A — 0. Thus we derive 
the shell crossing condition 



sinh 77 (sinh 77 — 77) 8 
(cosh 77 — l) 2 9 



(A8) 



Shell crossing first happens at 77,5c = 3.488 among the bound- 
ary shells, i.e., Ri — Ro in the above criterion. At shell cross- 
ing, the void interior has a relative density 



1 + A S 



9 (sinh 77^ - Vbc) 2 
2 (cosh77 sc -l) 3 



= 0.2047, 



(A9) 



which implies that the void has expanded by a factor of 
(1 + A sc ) -1 ' 3 = 1.697 in comoving radius. Note that these 
numbers do not depend on the size of the void. 

To calculate the linear theory prediction of the void 
underdensity at shell crossing <5 V , we expand R(t) to the 
first order with the help of the parametric solution (|A4|) 
and (|A"5)) , for R { < Ro 



R_ 
Ri 



a 
"i 



1 - 



So /3 



-Hit 



+ 



(A10) 



the first order of which gives the linear underdensity 



5 = So ( !#it 



•— [6(sinh 77-77)] 2/3 . 



(All) 



Thus the linear underdensity at shell crossing is 

5 v = S(r lsc ) = -2.717. (A12) 

Note this number is different in different cosmologies, e.g 
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Figure AI. Linear underdensity as a function of void expansion 
factor (exact: dot-dashed line, fit of equation (IB II ) : blue solid line). 
Also shown is the physicality constraint from requiring that the 
total volume fraction in voids from equation I I14H remain less than 
unity in SVdW model for the least stringent value of <5 C = 1.06 
(shaded region). Only as |<5 V | —> 0, where the regions are at the 
mean density and do not expand, does the model regain physi- 
cality. Bottom panel shows fractional error in the fit. 



it is S v — —2.731 for ACDM. However we find that such a 
small change in S v going from an EdS to a ACDM universe 
has a small impact on the predicted abundance of voids in 
the excursion set theor y. Also note the value for S v for the 
EdS universe quoted in ISheth fc van de Weygaert] (|2004h is 
—2.81 which was based on an approximate calculation and 
not on an exact treatment of the evolution of a tophat un- 
derdensitylj 

Similarly, for the spherical collapse model 



S = So I -Ht 



— [6(7/ -sin 77)] 2/3 . 



(A13) 



The well-known turn-around and virialization of halos occur 
at ?7ta = Ti" and ?7vir = 2-7T, leading to S c — 1.062 and S c — 
1.686 respectively. For the ACDM model these become 5 C = 
1.303 and S c = 1.674, for turnaround and collapse at z — 0. 
The EdS range encompasses that of the ACDM parameters 
and so in the main text we have adopted the EdS parameters 
to show the full range of possibilities. 

Finally note that these linear density thresholds 5 V and 
5 C , which are to be used in the excursion set formalism, are 
independent of the size of the structures. 



R. Sheth, private communication 
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APPENDIX B: SVDW MODEL 
MODIFICATIONS 

In this Appendix we explore modifications of the SVdW 
model prescription and appro ximations introduced in 
ISheth fc van de WevgaertJ (|2004 ). The spherical evolution 
model relates the linear underdensity S to nonlinear under- 
density A, or alternatively to r/rL = (1 -+- A)~3, where r is 
the void radius r — R(t; Ro) / a(t) and j"l the linear radius 
7~l = Ro/cii. If we relax the criterion for defining a void to 
correspond to underdense regions that have undergone shell 
crossing, there is additional freedom in defining the void 
abundance as a function of radius so long as 5 V and r/rL are 
chosen self-consistently. 

We show this relation before shell crossing for the EdS 
mode l in Fig. IA1I The relation is well fit to (|Bernardeaul 
1 19941 ) 



S v a c[l - (r/r L ) J/ 1 = c[l - (pv/p m ) H 



l/c 



(Bl) 



where c = 1.594, with errors below 0.2%. Also shown is the 
maximum r/rt,, constrained by requiring the total volume 
fraction in voids from equation (|14|) be less than 1, in the 
SVdW model. Clearly this constraint depends on S c , which 
in the plot is chosen to be the least stringent value 1.06 in 
the expected range. Note that no choice of <5 V and r/rL is 
physical for they all violate the total volume condition. In 
the main text, we also considered ad hoc modifications of 
the model where 5 V and r/rL a r e con sidered unrelated. 

ISheth fc van de WevgaertJ (|2004T ) also utilised an ap- 
proximation to the exact prediction for the abundance func- 
tion of equation ([6} which introduces notable errors for 
scales where the void-in-cloud process dominates and con- 
sequently the total volume fraction. Their approximation 

/w(cr )«y_ exp (__j ex p^____2_j, (B2) 

where v = S^/a 2 (M), had a stated realm of validity of 
Sc/\S V \ > 1/4 or T> < 4/5. Unfortunately, this approximation 
has uncontrolled errors at v <C 1, exactly where the void-in- 
cloud process operates as shown in Fig. IB1I Our piecewise 
approximation is accurate at the 0.2% level or better every- 
where. Note that the errors and smoothness of our approxi- 
mation can be improved at the transition point by a suitable 
interpolation between the two piecewise curves though it is 
not necessary for this work. 
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